Clustering and Classification of Satellite Imagery

Abstract

-RAQUEL # Unsupervised Clustering -RAQUEL any intro stuff about this topic, feell free to use references This section will use unsupervised clustering to split up given satellite images into several groups, with the goal of identifying meaningfull clusters that represent different geographic features such as water, farm land, forest or buildings.

## Loading required package: raster
## Loading required package: sp
## Loading required package: rgdal
## rgdal: version: 1.1-1, (SVN revision 572)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
##  Path to GDAL shared files: C:/Program Files/R/R-3.2.3/library/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: C:/Program Files/R/R-3.2.3/library/rgdal/proj
##  Linking to sp version: 1.2-1 
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
## [1] 24000
#rasterJ<-brick("C:\\Users\\Owner\\Documents\\Portugal\\Sat_Image_Classification\\data\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
rasterJ<-brick("C:\\Users\\D059331\\Desktop\\DM GIC\\data\\img\\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF")
rasterJ
## class       : RasterBrick 
## dimensions  : 4238, 4239, 17964882, 4  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -105092.3, -94918.74, -50088.84, -39917.64  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : C:\Users\D059331\Desktop\DM GIC\data\img\J_04AUG14112729-M2AS-000000137917_01_P001_etrs89.TIF 
## names       : J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3, J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4 
## min values  :                                                  0,                                                  0,                                                  0,                                                  0 
## max values  :                                                255,                                                255,                                                255,                                                255
plotRGB(rasterJ, 3,2,1)

ext <- extent(-104294.4, -102964.5, -43623.48, -42742.44 )
sector <- crop(rasterJ, ext)
plotRGB(sector, 3, 2, 1)

sector[[5]] <- normdiff(sector)
plot(sector)

summary(sector)
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.1
## Min.                                                    23
## 1st Qu.                                                 31
## Median                                                  34
## 3rd Qu.                                                 40
## Max.                                                   176
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.2
## Min.                                                    26
## 1st Qu.                                                 40
## Median                                                  50
## 3rd Qu.                                                 65
## Max.                                                   255
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.3
## Min.                                                     8
## 1st Qu.                                                 17
## Median                                                  34
## 3rd Qu.                                                 51
## Max.                                                   194
## NA's                                                     0
##         J_04AUG14112729.M2AS.000000137917_01_P001_etrs89.4      values
## Min.                                                     4 -0.64885496
## 1st Qu.                                                 10 -0.20000000
## Median                                                  53  0.05645161
## 3rd Qu.                                                 70  0.15384615
## Max.                                                   191  0.76870748
## NA's                                                     0  0.00000000
means<-mov.fun(sector,3,mean, normalize=FALSE, verbose=TRUE)
## >> Band  1  of  5  done.
## >> Band  2  of  5  done.
## >> Band  3  of  5  done.
## >> Band  4  of  5  done.
## >> Band  5  of  5  done.
sds<-mov.fun(sector,3,sd, normalize=FALSE, verbose=TRUE)
## >> Band  1  of  5  done.
## >> Band  2  of  5  done.
## >> Band  3  of  5  done.
## >> Band  4  of  5  done.
## >> Band  5  of  5  done.
stk<-stack(means,sds)
sectorMovings<-brick(stk)
summary(sectorMovings)
##         layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1       layer.1
## Min.       24.00000    27.88889    9.111111    4.555556   -0.51877395
## 1st Qu.    30.88889    40.44444   16.666667   10.111111   -0.21001302
## Median     34.55556    51.00000   35.555556   53.555556    0.05963125
## 3rd Qu.    40.22222    64.88889   50.888889   70.222222    0.16173571
## Max.      172.11111   255.00000  189.222222  171.111111    0.73729629
## NA's     1840.00000  1840.00000 1840.000000 1840.000000 1840.00000000
##          layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2      layer.2
## Min.       0.0000000    0.000000    0.000000    0.000000 0.000000e+00
## 1st Qu.    0.9718253    1.666667    1.322876    1.414214 2.197058e-02
## Median     1.4142136    2.877113    2.438123    2.915476 4.184848e-02
## 3rd Qu.    2.8037673    6.009252    5.967505    7.270565 6.507355e-02
## Max.      53.2846132   85.835275   64.031242   60.055761 3.334456e-01
## NA's    1840.0000000 1840.000000 1840.000000 1840.000000 1.840000e+03
sectorMovings
## class       : RasterBrick 
## dimensions  : 367, 555, 203685, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -104295.5, -102963.5, -43623.24, -42742.44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1,     layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2,     layer.2 
## min values  :  24.0000000,  27.8888889,   9.1111111,   4.5555556,  -0.5187739,   0.0000000,   0.0000000,   0.0000000,   0.0000000,   0.0000000 
## max values  : 172.1111111, 255.0000000, 189.2222222, 171.1111111,   0.7372963,  53.2846132,  85.8352751,  64.0312424,  60.0557611,   0.3334456
plot(sectorMovings)

sectorMovingsTrim<-trim(sectorMovings, padding=0, values=NA)
summary(sectorMovingsTrim)
##         layer.1.1.1 layer.2.1.1 layer.1.2.1 layer.2.2.1     layer.1
## Min.       24.00000    27.88889    9.111111    4.555556 -0.51877395
## 1st Qu.    30.88889    40.44444   16.666667   10.111111 -0.21001302
## Median     34.55556    51.00000   35.555556   53.555556  0.05963125
## 3rd Qu.    40.22222    64.88889   50.888889   70.222222  0.16173571
## Max.      172.11111   255.00000  189.222222  171.111111  0.73729629
## NA's        0.00000     0.00000    0.000000    0.000000  0.00000000
##         layer.1.1.2 layer.2.1.2 layer.1.2.2 layer.2.2.2    layer.2
## Min.      0.0000000    0.000000    0.000000    0.000000 0.00000000
## 1st Qu.   0.9718253    1.666667    1.322876    1.414214 0.02197058
## Median    1.4142136    2.877113    2.438123    2.915476 0.04184848
## 3rd Qu.   2.8037673    6.009252    5.967505    7.270565 0.06507355
## Max.     53.2846132   85.835275   64.031242   60.055761 0.33344559
## NA's      0.0000000    0.000000    0.000000    0.000000 0.00000000
plot(sectorMovingsTrim)

sectorMovingsTrimNorm<-raster.scale(sectorMovingsTrim)
sectorMovingsTrimNorm
## class       : RasterBrick 
## dimensions  : 365, 553, 201845, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 2.4, 2.4  (x, y)
## extent      : -104293.1, -102965.9, -43620.84, -42744.84  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=39.66825833333333 +lon_0=-8.133108333333334 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer.1.1.1, layer.2.1.1, layer.1.2.1, layer.2.2.1, layer.1, layer.1.1.2, layer.2.1.2, layer.1.2.2, layer.2.2.2, layer.2 
## min values  :           0,           0,           0,           0,       0,           0,           0,           0,           0,       0 
## max values  :           1,           1,           1,           1,       1,           1,           1,           1,           1,       1
##           used (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 1043532 55.8    1770749  94.6  1770749  94.6
## Vcells 6149593 47.0   18416227 140.6 20881747 159.4
kmeansCenters <- performKMeans(sectorMovingsTrimNorm, 50, output="centers")
interClusterDistance <- dist(kmeansCenters)
clusterHierarchy <- hclust(interClusterDistance, method="ward.D")
plot(clusterHierarchy)

- decide to use 8 clusters (graphical argument)

clusterRaster <- performKMeans(sectorMovingsTrimNorm, 8)
plot(clusterRaster)